# Replication materials for 
# Peisker, Hoffmann, Muttarak (2025)
# Climate news mediates extreme weather effects on climate change concern

packages <- c(
  "tidyverse", "data.table", "countrycode", "scales",
  "fixest", "parallel", "rnaturalearth", "Cairo",
  "texreg", "xtable", "sf", "readxl",
  "scico", "patchwork"
)
lapply(packages, library, character.only = TRUE)
lapply(packages, citation)

source("R Files/functions.R")
theme_set(theme_bw())
setDTthreads(parallel::detectCores())

# #### load data ####
load("Data/eb_climate_weekly.RData")
load("Data/final eurobarometer monthly.RData")

#### plot theme ####
theme_set(theme_bw())
map_pal <- "vik"
map_dir <- -1
map_begin <- 0
map_end <- 0.5

trend_pal <- "vik"
trend_dir <- 1
trend_begin <- 0.2
trend_end <- 0.8

unique(eb_climate_weekly$nuts0)
length(unique(eb_climate_weekly$nuts0))
length(unique(eb_climate_weekly$nuts))
length(unique(paste(eb_climate_weekly$year, eb_climate_weekly$month)))

regions_raw <- st_read("Data/NUTS_RG_10M_2016_3035.shp/NUTS_RG_10M_2016_3035.shp")
countries_raw <- 
  regions_raw %>% filter(LEVL_CODE == 0)
other_countries <- ne_countries(
  scale = 50, 
  country = c(
    "Tunisia", "Algeria", "Morocco", "Turkey",
    "Bosnia and Herzegovina", "Moldova", "Kosovo",
    "Moldova", "Ukraine", "Russia", "Belarus",
    "Georgia", "Syria","Iraq", "Azerbaijan", "Iran", "Lebanon",
    "Israel", "Saudi Arabia", "Palestine", "Jordan",
    "Armenia", "Kazakhstan", "Uzbekistan", "Turkmenistan"
  ),
  returnclass = "sf"
)

#### media means over time and bar chart ####
emm_vars <- c("emm_climate_w_nat_4w", "emm_climate_w_reg_4w")
emm_rmean <- 
  eb_climate_weekly %>% 
  melt(id.vars = c("t", "year", "week", "nuts", "nuts0", "region"), measure.vars = emm_vars) %>% 
  .[,`:=`( 
    Region = region,
    Topic = factor(fcase(
      grepl("climate", variable), "Climate change",
      grepl("env", variable), "Environment & ecology"
      ), levels = c("Climate change", "Environment & ecology")
    ),
    Publication = factor(case_when(
      grepl("_tot", variable, fixed = TRUE) ~ "All",
      grepl("_reg", variable, fixed = TRUE) ~ "Regional",
      grepl("_nat", variable, fixed = TRUE) ~ "National"
    ), levels = c("All", "Regional", "National")),
    t = year + week / 53
  )] %>% 
  .[,
    .(emm_mean = mean(value)),
    by = .(t, Topic, Publication)#
  ]

events_plot <- 
  read_csv("Data/events_eu.csv") %>% 
  filter(# filter some overlapping events for readability
    plot == 1
  ) %>% 
  mutate(
    date_start = as.Date(paste(year_start, month_start, day_start, sep = "-")) - 4,
    week_start = week(date_start),
    date_end = as.Date(paste(year_end, month_end, day_end, sep = "-")) + 4,
    week_end = week(date_end),
    start = as.integer(year_start) + as.integer(week_start) / 53,
    end = as.integer(year_end) + as.integer(week_end) / 53,
    middle = (start + end) / 2
  ) %>%
  group_by(event_name) %>% 
  summarise(across(c("start", "end", "middle"), unique)) %>%
  arrange(middle) %>% ungroup %>% 
  mutate(
    event_number = 1:nrow(.)
  )

# for figure caption
paste0("(", events_plot$event_number, ") ", events_plot$event_name, collapse = ", ")

emm_trends <- 
  emm_rmean[Topic == "Climate change" & Publication %in% c("Regional", "National"), ] %>% 
  ggplot(aes(x = t, y = emm_mean)) +
  geom_rect(
    aes(xmin = start, xmax = end, ymin = 0, ymax = max(emm_rmean$emm_mean, na.rm = TRUE), group = event_name),
    data = events_plot, inherit.aes = FALSE, alpha = 0.3
  ) +
  geom_line(aes(color = Publication), linewidth = 0.8) +
  geom_label(
    data = events_plot, aes(
      x = middle,
      y = c(rep(c(max(emm_rmean$emm_mean, na.rm = TRUE) - 0.3, max(emm_rmean$emm_mean, na.rm = TRUE) - 0.6), length.out = nrow(events_plot))),
      label = event_number)
  ) +
  labs(y = "Climate news (4-week mean)", color = "Publications") +
  scale_x_continuous(limits = c(2014.4, 2020.1), breaks = 2015:2020, expand = c(0, 0)) +#
  scale_y_continuous(expand = c(0, 0)) +
  theme(
    legend.position = "bottom",
    legend.box="vertical", #legend.margin=margin()
    axis.title.x = element_blank()
  ) +
  scale_color_scico_d(
    palette = trend_pal, direction = trend_dir,
    begin = 0.2, end = 0.8
  ) + 
  theme(text = element_text(size = 18))
emm_trends

#### scatter plot ####
emm_scatter <- 
  eb_climate_weekly %>% 
  filter(!is.na(cc_eu_concern_w)) %>% 
  pivot_longer(c(emm_climate_w_reg_4w, emm_climate_w_nat_4w), names_to = "level", values_to = "news") %>% 
  mutate(
    publication = factor(
      if_else(grepl("_reg", level), "Regional", "National"),
      levels = c("Regional", "National")
      )
  ) %>% 
  ggplot(aes(x = cc_eu_concern_w, y = news, color = publication)) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm", se = F) +
  scale_x_continuous(labels = label_percent()) +
  labs(x = "Climate concern", y = "Climate news (4-week mean)", color = "Publications") +
  scale_color_scico_d(
    palette = trend_pal, direction = trend_dir,
    begin = 0.2, end = 0.8
  ) + 
  theme(text = element_text(size = 18), legend.position = "bottom")
emm_scatter

#### concern over time ####
eb_plot <- 
  eb_month %>% 
  mutate(
    nuts0 = substr(nuts, 1, 2),
    Region = case_when(
      nuts0 %in% c("DK", "FI", "IS", "SE", "UK", "IE") ~ "North",
      nuts0 %in% c("AT", "BE", "DE", "FR", "LU", "NL") ~ "West",
      nuts0 %in% c("EL", "CY", "ES", "IT", "MT", "PT") ~ "South",
      nuts0 %in% c("AL", "BG", "CZ", "HU", "HR", "EE", "LT", "LV", "ME", "MK", "PL", "RO", "RS", "SI", "SK", "TR") ~ "East"
    )) %>% 
  group_by(year, Region) %>% 
  summarise(
    env_concern_w = mean(env_concern_w, na.rm = TRUE) * 100,
    cc_eu_concern_w = mean(cc_eu_concern_w, na.rm = TRUE) * 100,
    cc_world_concern_w = mean(cc_world_concern_w, na.rm = TRUE) * 100
  )

concern_line <- 
  eb_plot %>% 
  ggplot(aes(x = year, y = cc_eu_concern_w)) +
  annotate("rect", xmin = 2014, xmax = 2019, ymin = 0, ymax = 34, alpha = 0.2) +
  geom_line(aes(color = Region), linewidth = 0.75) +
  geom_point(aes(shape = Region, color = Region), size = 2) +
  scale_x_continuous(limits = c(2010,2019), breaks = seq(2000, 2020, 1)) +
  scale_y_continuous(
    limits = c(-0, 34), breaks = seq(0, 30, 5), expand = c(0, 0),
    labels = label_percent(scale = 1)
  ) +
  scale_color_scico_d(
    palette = trend_pal, direction = trend_dir,
    begin = trend_begin, end = trend_end
  ) + 
  labs(x = "Year", y = "Climate change concern") +
  theme(
    text = element_text(size = 18),
    legend.position = "bottom",
    axis.title.x = element_blank()
  )
concern_line

eb_plot_map <- 
  eb_month %>% 
  filter(between(year, 2014, 2019)) %>% 
  group_by(nuts) %>% 
  summarise(
    env_concern_w = mean(env_concern_w, na.rm = TRUE) * 100,
    cc_eu_concern_w = mean(cc_eu_concern_w, na.rm = TRUE) * 100,
    cc_world_concern_w = mean(cc_world_concern_w, na.rm = TRUE) * 100
  )

concern_map <- 
  regions_raw %>% 
  merge(eb_plot_map, by.y = "nuts", by.x = "NUTS_ID") %>% 
  ggplot() +
  geom_sf(data = countries_raw, lwd = 0, fill = "grey80") +
  geom_sf(data = other_countries, size = 0, fill = "grey80") +
  geom_sf(aes(fill = cc_eu_concern_w), colour = NA) +
  geom_sf(data = countries_raw, lwd = 0.3, fill = NA) +
  geom_sf(data = other_countries, size = 0.3, fill = NA) +
  scale_fill_scico(
    name = "Climate change concern \n(mean 2014-2019) ",
    labels = label_percent(scale = 1),
    breaks = seq(0, 100, 10),
    na.value = NA,
    palette = map_pal, direction = map_dir,
    begin = map_begin, end = map_end
  ) + #
  coord_sf(
    crs = 3035, #4326
    xlim = c(240e4, 630e4),#-12, 34
    ylim = c(130e4, 550e4),#33, 73, 
    expand = FALSE
  ) +
  theme(text = element_text(size = 18), legend.position = "bottom")
concern_map

emm_trends / (concern_line + emm_scatter) +
  plot_annotation(tag_levels = "a") 
ggsave(filename = "Figures/media_concern_trend.pdf", width = 14, height = 12)

#### google trends ####
google_plot <- 
  eb_climate_weekly %>% 
  # filter UK for trend because it is an outlier 
  filter(!(nuts %in% c("UKC", "UKD", "UKE", "UKF", "UKG", "UKH", "UKI", "UKJ", "UKK"))) %>% 
  .[,
    .(
      google_mean = mean(google_trend_4w, na.rm = TRUE),
      time = year + week / 53
      ),
    by = .(region, year, week)#
  ]

google_line <- 
  google_plot %>% 
  ggplot(aes(x = time, y = google_mean)) +
  geom_line(aes(color = region), linewidth = 0.75) +
  scale_x_continuous(breaks = seq(2000, 2020, 1)) +
  scale_y_continuous(
    limits = c(0, NA), breaks = seq(0, 30, 5)#, expand = c(0, 2)
  ) +
  scale_color_scico_d(
    palette = trend_pal, direction = trend_dir,
    begin = trend_begin, end = trend_end
  ) + 
  labs(
    x = "Year", 
    y = "Google climate searches (4-week mean)", 
    color = "Region"
    ) +
  theme(
    legend.position = "bottom",
    axis.title.x = element_blank()
  )
google_line

google_plot_map <- 
  eb_climate_weekly %>% 
  filter(between(year, 2014, 2019)) %>% 
  group_by(nuts) %>% 
  summarise(
    google_trend = mean(google_trend, na.rm = TRUE)
  )

google_map <- 
  regions_raw %>% 
  merge(google_plot_map, by.y = "nuts", by.x = "NUTS_ID") %>% 
  ggplot() +
  geom_sf(data = countries_raw, lwd = 0, fill = "grey80") +
  geom_sf(data = other_countries, size = 0, fill = "grey80") +
  geom_sf(aes(fill = google_trend), colour = NA) +
  geom_sf(data = countries_raw, lwd = 0.3, fill = NA) +
  geom_sf(data = other_countries, size = 0.3, fill = NA) +
  scale_fill_scico(
    name = "Google climate searches \n(regional mean)",
    limits = c(0, NA),
    breaks = seq(0, 100, 10),
    na.value = NA,
    palette = map_pal, direction = map_dir,
    begin = map_begin, end = map_end
  ) + #
  coord_sf(
    crs = 3035, #4326
    xlim = c(240e4, 630e4),#-12, 34
    ylim = c(130e4, 550e4),#33, 73, 
    expand = FALSE
  ) +
  theme(legend.position = "bottom")
google_map

google_line + google_map +
  plot_annotation(tag_levels = "a")
ggsave(filename = "Figures/google_trend.pdf", width = 9, height = 4.5)

#### summary table media ####
tab_weekly <-  read_csv("Tables/var_table_weekly.csv") 
in_vars <- 
  tab_weekly %>% 
  filter(!is.na(var)) %>% 
  select(var) %>% 
  .$var
round_cols <- c("Min", "25\\%", "50\\%", "75\\%", "Max", "SD", "Within SD")

within_sd <- 
  sapply(in_vars, function(i){ 
    felm_tmp <- feols(  
      f_feols(yvar = "tmean_nf", xvars = i, fevars = m_fv), 
      eb_climate_weekly[!is.na(emm_climate_w_reg_l1),], demeaned = TRUE)
    sd_tmp <- c(sd(felm_tmp$X_demeaned, na.rm = TRUE))
    return(sd_tmp)# %>% round(2)
  }) %>% 
  data.frame() %>% 
  rownames_to_column
names(within_sd) <- c("var", "Within SD")
within_sd

sum_table_media <- 
  eb_climate_weekly[!is.na(emm_climate_w_reg_l1), ..in_vars] %>% 
  melt(variable.name = "var") %>% 
  .[
    , .(
      N = sum(!is.na(value)),
      Min = min(value, na.rm = TRUE),
      `25\\%` = quantile(value, probs = 0.25, na.rm = TRUE),
      `50\\%` = quantile(value, probs = 0.50, na.rm = TRUE),#mean(value, na.rm = TRUE)
      `75\\%` = quantile(value, probs = 0.75, na.rm = TRUE),
      Max = max(value, na.rm = TRUE),
      SD = sd(value, na.rm = TRUE)
    ),
    by = .(var)
  ] %>% 
  merge(within_sd, by = "var", all.x = TRUE) %>% 
  .[
    , c(round_cols) := lapply(.SD, function(x){round(x, 3)}),
    .SDcols = c(round_cols)
  ]
print(sum_table_media)

tab_weekly2 <- 
  tab_weekly %>% 
  left_join(sum_table_media, by = "var")
tab_weekly2[ , c("Variable", "Source", "N", round_cols)] %>% print(n = 50)

print(
  xtable(
    tab_weekly2[ , c("Variable", "Source", "N", round_cols)],
    align = c(  "l",   "l",       "l",     "r", "r",    "r",    "r",      "r",    "r",   "r",   "r")
  ), #
  file = "Tables/summary_weekly.tex",
  include.rownames = FALSE,
  booktabs = TRUE,
  floating = FALSE,
  math.style.negative = TRUE,
  sanitize.text.function	= identity
)

#### summary table concern ####
tab_eb <-  read_csv("Tables/var_table_eb.csv") 
in_vars <- 
  tab_eb %>% 
  filter(!is.na(var)) %>% 
  select(var) %>% 
  .$var
round_cols <- c("Min", "25\\%", "50\\%", "75\\%", "Max", "SD", "Within SD")

within_sd_eb <- 
  sapply(in_vars, function(i){ 
    felm_tmp <- feols(  
      f_feols(yvar = "tmean_nf", xvars = i, fevars = m_fv), 
      eb_climate_weekly[!is.na(cc_eu_concern_w),], demeaned = TRUE)
    sd_tmp <- c(sd(felm_tmp$X_demeaned, na.rm = TRUE))
    return(sd_tmp)# %>% round(2)
  }) %>% 
  data.frame() %>% 
  rownames_to_column
names(within_sd_eb) <- c("var", "Within SD")
within_sd_eb

sum_table_concern <- 
  eb_climate_weekly[!is.na(cc_eu_concern_w), ..in_vars] %>% 
  melt(variable.name = "var") %>% 
  .[
    , .(
      N = sum(!is.na(value)),
      Min = min(value, na.rm = TRUE),
      `25\\%` = quantile(value, probs = 0.25, na.rm = TRUE),
      `50\\%` = quantile(value, probs = 0.50, na.rm = TRUE),#mean(value, na.rm = TRUE)
      `75\\%` = quantile(value, probs = 0.75, na.rm = TRUE),
      Max = max(value, na.rm = TRUE),
      SD = sd(value, na.rm = TRUE)
    ),
    by = .(var)
  ] %>% 
  merge(within_sd_eb, by = "var", all.x = TRUE) %>% 
  .[
    , c(round_cols) := lapply(.SD, function(x){round(x, 3)}),
    .SDcols = c(round_cols)
  ]
print(sum_table_concern)

tab_eb2 <- 
  tab_eb %>% 
  left_join(sum_table_concern, by = "var")
tab_eb2[ , c("Variable", "Source", "N", round_cols)] %>% print(n = 50)

print(
  xtable(
    tab_eb2[ , c("Variable", "Source", "N", round_cols)],
    align = c(  "l",   "l",       "l",     "r", "r",    "r",    "r",      "r",    "r",   "r",   "r")
  ), #
  file = "Tables/summary_eb.tex",
  include.rownames = FALSE,
  booktabs = TRUE,
  floating = FALSE,
  math.style.negative = TRUE,
  sanitize.text.function = identity
)

#### data summary by country ####
emm_sources <- read_csv("Data/tab_emm_sources.csv", col_types = c("c","i","i"))
emm_sources

cc_concern <- 
  eb_climate_weekly[
    !is.na(cc_eu_concern_w), 
    .(
      Country = unique(countrycode(nuts0, "eurostat", "country.name")),
      # N = .N,
      `Macro-region` = unique(region),
      `NUTS level` = min(fcase(
        nchar(nuts) == 2, 0L,
        nchar(nuts) == 3, 1L,
        nchar(nuts) == 4, 2L,
        nchar(nuts) == 5, 3L
      )),
      `NUTS regions` = length(unique(nuts)),
      `Eurobarometer waves` = as.integer(max(table(nuts)))
    ),
    by = .(nuts0)
  ] %>% 
  merge(emm_sources) %>% 
  .[order(Country), -2] 

print(
  xtable::xtable(cc_concern, digits = 0), 
  file = "Tables/data_summary.tex",
  # include.rownames = FALSE,
  booktabs = TRUE,
  floating = FALSE,
  math.style.negative = TRUE,
  sanitize.text.function	= identity
)

#### list of events ####
events_table <- 
  read_csv("Data/events_eu.csv") %>% 
  mutate(
    date_start = as.Date(paste(year_start, month_start, day_start, sep = "-")),
    date_end = as.Date(paste(year_start, month_end, day_end, sep = "-"))
  ) %>% 
  filter(
    date_start > "2014-05-23",
    date_end < "2019-12-31"
  ) %>%
  arrange(date_start) %>% 
  transmute(
    `Event name` = event_name,
    `Start date` = paste(day_start, month_start, year_start, sep = "."),
    `End date` = paste(day_end, month_end, year_end, sep = ".")
  )
events_table %>% print(n = 50)
print(
  xtable(events_table, align = c("l", "l", "r", "r")), #
  file = "Tables/tab_events_list.tex",
  include.rownames = TRUE,
  booktabs = TRUE,
  floating = FALSE,
  math.style.negative = TRUE,
  sanitize.text.function	= identity
)

